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ABSTRACT 

We study the dynamical evolution of globular clusters using our 2D Monte Carlo code with the 
inclusion of primordial binary interactions for equal-mass stars. We use approximate analytical cross 
sections for energy generation from binary-binary and binary-single interactions. After a brief period 
of slight contraction or expansion of the core over the first few relaxation times, all clusters enter a 
much longer phase of stable "binary burning" lasting many tens of relaxation times. The structural 
parameters of our models during this phase match well those of most observed globular clusters. At 
the end of this phase, clusters that have survived tidal disruption undergo deep core collapse, followed 
by gravothermal oscillations. Our results clearly show that the presence of even a small fraction of 
binaries in a cluster is sufficient to support the core against collapse significantly beyond the normal 
core collapse time predicted without the presence of binaries. For tidally truncated systems, collapse is 
easily delayed sufficiently that the cluster will undergo complete tidal disruption before core collapse. 
As a first step toward the eventual goal of computing all interactions exactly using dynamical three- 
and four-body integration, we have incorporated an exact treatment of binary-single interactions in our 
code. We show that results using analytical cross sections are in good agreement with those using exact 
three-body integration, even for small binary fractions where binary-single interactions arc energetically 
most important. 

Subject headings: celestial mechanics, stellar dynamics — globular clusters: general — methods: 
numerical 



1. INTRODUCTION 

The realization about 10 years ago that primordial bi- 
naries are present in globular clusters in dynamically sig- 
nificant numbers has completely changed our theoretical 
perspective on these systems (see, e.g., the review by Hut 
et al. 1992a). Most importantly, dynamical interactions 
between hard primordial binaries and other single stars or 
binaries are thought to be the primary energy generation 
mechanism responsible for supporting a globular cluster 
against core collapse (Goodman & Hut 1989; McMillan 
et al. 1990, 1991; Gao et al. 1991). The term "binary burn- 
ing" is now often used by analogy with hydrogen burning 
for stars. In the same way that hydrogen burning allows a 
star like the Sun to remain in thermal equilibrium on the 
main sequence for a time much longer than the Kelvin- 
Hclmholtz timescale, primordial binary burning allows a 
globular cluster to maintain itself in quasi-thermal equilib- 
rium and avoid core collapse for a time much longer than 
the two-body relaxation timescale. 

In addition, strong dynamical interactions involving bi- 
naries can explain very naturally the large numbers of ex- 
otic objects found in dense star clusters. Exchange in- 
teractions between hard primordial binaries and neutron 
stars inevitably produce large numbers of X-ray binaries 
and recycled pulsars in globular clusters (Hut et al. 1991; 
Sigurdsson & Phinney 1995; Davies & Hansen 1998; Rasio 
et al. 2000). Resonant interactions of primordial binaries 
result in dramatically increased collision rates for main- 
sequence stars in globular clusters and even open clus- 



ters (Bacon et al. 1996; Cheung et al. 2003; Leonard 1989; 
Leonard & Linnell 1992). Direct observational evidence 
for stellar collisions and mergers of main-sequence stars in 
globular clusters comes from the detection of large num- 
bers of bright blue stragglers concentrated in the dense 
cluster cores (Bailyn 1995; Bellazzini et al. 2002; Ferraro 
et al. 2001). Previously it was thought that primordial 
binaries were essentially nonexistent in globular clusters, 
and so other mechanisms such as tidal capture and three- 
body encounters had to be invoked in order to form bi- 
naries dynamically during deep core collapse. However, 
these other mechanisms have some serious problems, and 
are much more likely to result in mergers than in the 
formation of long-lived binaries (Chernoff & Huang 1996; 
Kochanek 1992; Kumar & Goodman 1996; Portegies Zwart 
& McMillan 2002; Kim & Lee 1999; Kim et al. 1998; Lee & 
Ostriker 1993). Multiple mergers of main-sequence stars 
and runaway collisions in young star clusters could lead 
to the formation of a massive central black hole in some 
systems (Lee 1993; Gebhardt et al. 2002; Portegies Zwart 
& McMillan 2002). 

The primordial binary fraction is therefore a key input 
parameter for any realistic study of dense star cluster dy- 
namics (Hut et al. 1992a). Early determinations of bi- 
nary fractions in globular clusters came from observations 
of spectroscopic binaries with red giant primaries (Pryor 
et al. 1988, see, e.g., Cote et al. 1996 for a more recent 
study) as well as eclipsing binaries (Mateo et al. 1990; Yan 
& Mateo 1994). Hubble Space Telescope observations have 
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provided direct constraints on primordial binary fractions 
in the central regions of many globular clusters, where bi- 
naries are expected to concentrate because of mass segre- 
gation. Rubenstein & Bailyn (1997) used observations of a 
broadened main sequence in NGC 6752 to derive a binary 
fraction in the range 15%-38% for the inner cluster core. 
Their method has now been applied to many other clus- 
ters. For example, Bellazzini et al. (2002) derive a similar 
binary fraction, in the range 0.08-0.38, in the central re- 
gion of NGC 288. Adding proper motion information can 
lead to much tighter constraints, as in the case of NGC 
6397, where Cool & Bolton (2002) derive a binary fraction 
< 5 — 7% near the center. 

Despite the crucial role of primordial binaries in the dy- 
namical evolution of a dense star cluster, the overall evo- 
lution of the binary population within a cluster, and its 
direct implications for the formation rate of observable 
systems such as recycled pulsars and blue stragglers, re- 
mains poorly understood theoretically. One reason is that 
the relative importance of binary interactions in a cluster, 
like many other dynamical processes, depends in a complex 
manner on the number of stars in the system. This makes 
it difficult to extend results obtained from small direct A- 
body simulations to realistic globular cluster models. In 
particular, the rate at which binaries are "burned" and, 
ultimately, destroyed or ejected from the cluster depends 
on the size of the cluster. When the initial primordial bi- 
nary fraction is below a certain critical value, a globular 
cluster core can run out of binaries before the end of its life- 
time, i.e., before disruption in the tidal field of the Galaxy 
(McMillan & Hut 1994). Without the support of binaries, 
the cluster will then undergo a much deeper core collapse, 
perhaps followed by gravothermal oscillations (Sugimoto 
& Bettwieser 1983; Breeden et al. 1994; Makino 1996). At 
maximum contraction, the core density may increase by 
many orders of magnitude, leading to greatly enhanced 
interaction rates. 

Detailed numerical studies of globular cluster evolution 
with primordial binaries are still lacking, for several rea- 
sons. First, the inclusion of even a modest fraction of 
primordial binaries adds a very significant computational 
overhead to A-body simulations. This is mainly due to 
the extra computations required to treat binary interac- 
tions, but also because the lifetime of a cluster can be 
significantly extended (by up to many orders of magni- 
tude) through binary burning. In addition, in direct N- 
body simulations, the extremely large ratio of the overall 
cluster dynamical time to the orbital period of close bina- 
ries (as large as ~ 10 10 in a globular cluster!) introduces 
many computational difficulties. This makes A-body sim- 
ulations with primordial binaries prohibitively expensive 
for A > 10 4 stars, although special-purpose supercom- 
puters such as the new GRAPE-6 may increase this limit 
in the near future (Makino 2001). Orbit-averaged calcu- 
lations, like direct Fokker-Planck integrations and Monte 
Carlo simulations, get around this problem by treating bi- 
naries just like single stars, except during brief periods 
of strong interactions. Unfortunately, this requires that 
cross sections for strong interactions involving binaries be 
known accurately, for a wide range of binary parameters 
(masses, semi-major axes, and eccentricities). These cross 
sections are difficult to determine in general, and reliable 



semi-analytic fits to numerical scattering experiments are 
only available for simple configurations such as those in- 
volving equal-mass stars. For these reasons, most previ- 
ous numerical studies of globular clusters with primordial 
binaries have been limited either to clusters with equal- 
mass stars (Gao et al. 1991; Spitzer & Mathieu 1980), or 
to very small clusters with A ~ 10 3 — 10 4 stars (Heggie 
& Aarseth 1992; Hurley et al. 2001; McMillan et al. 1990, 
1991; McMillan & Hut 1994). Simplified treatments have 
also been employed in which the dynamics of the bina- 
ries was followed in a static cluster background (Hut et al. 
1992b) or in a background cluster modeled as an evolving 
gas sphere (Giersz & Spurzem 2000). 

The results of Gao et al. (1991, hereafter GGCM91), 
based on direct Fokker-Planck integrations, were the first 
to clearly illustrate the dominant effect of even a small 
fraction of primordial binaries on the evolution of a glob- 
ular cluster. In this paper, we present the first study of 
globular cluster evolution with primordial binaries based 
on self-consistent Monte Carlo simulations with a realis- 
tically large number of stars (A > 10 5 ). Partly in or- 
der to allow better comparison of our results with those 
of GGCM91, we use similar initial conditions and cross 
sections for binary-binary and binary-single interactions, 
even though our method for implementing these cross sec- 
tions in the Monte Carlo scheme is completely different. 
In addition, the results of GGCM91 were obtained using 
a 1-D Fokker-Planck method (in which isotropy in veloc- 
ity space is enforced). More realistic 2-D (anisotropic) 
Fokker-Planck calculations with primordial binaries have 
never been reported in the literature, to the best of our 
knowledge. Even for the 1-D calculations, and with only 
a single parameter representing the internal structure of 
binaries (namely, their binding energy), the inclusion of 
binary-binary interactions significantly increased the over- 
all computation time. Since the Fokker-Planck method 
uses distribution functions to represent the system, every 
new parameter adds a new dimension to the phase space, 
making the Fokker-Planck equation more difficult to solve 
numerically. It has also been shown recently that the 1-D 
treatment is inadequate in dealing with some aspects of 
the evolution, such as the escape rate from tidally trun- 
cated clusters (Takahashi & Portegies Zwart 1998, 2000). 
Many difficulties in the direct Fokker-Planck approach 
come from the basic representation of the system in terms 
of smooth distribution functions. Neglecting the discrete 
nature of the system makes it impossible to follow the de- 
tails of individual interactions, such as binary-single or 
binary-binary interactions. The implicit assumption that 
A — > oo also makes it difficult to scale the results of direct 
Fokker-Planck simulations to finite systems with different 
numbers of stars. 

Our Monte Carlo method provides an intermediate ap- 
proach, which combines many of the benefits of direct A- 
body simulations (such as the description of the cluster on 
a star-by-star basis and the possibility to treat each in- 
dividual interaction in detail) with the speed of an orbit- 
averaged calculation. Our method is also 2-D in velocity 
space by construction, and hence properly accounts for 
any velocity anisotropy in the system. Another benefit of 
the method is that it allows a wide range of binary param- 
eters to be used without having to modify the underlying 
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orbit-averaged calculation of the relaxation processes. In 
principle, individual interactions can be treated in as much 
detail as in direct iV-body simulations, by computing all 
strong encounters exactly using three-body or four-body 
integrators. As a first step, for this paper, we have incor- 
porated a three-body dynamical integrator into our code, 
which allows binary-single interactions to be computed 
exactly (without reference to approximate, pre-compilcd 
cross sections). This allows us to follow the outcomes of 
interactions more precisely, and, most importantly, will 
allow us in the future to extend our code to multi-mass 
systems, for which analytic cross sections are not avail- 
able. 

2. TREATMENT OF BINARY INTERACTIONS 

We use the basic Henon-type Monte Carlo method for 
modeling the dynamical evolution of clusters as a sequence 
of equilibrium models subject to regular velocity pertur- 
bations (Henon 1971a, b); our code has been described in 
detail by Joshi et al. (2000, 2001, hereafter Papers I and 
II). The regular velocity perturbations are calculated us- 
ing Hcnon's method to represent the average effect of many 
long-range small-angle gravitational scattering encounters 
using one suitably chosen encounter with a nearby star 
(Henon 1971b). At each time step, we calculate the Monte- 
Carlo realized radial position and velocity of each star (as- 
suming spherical symmetry), which we use to calculate 
whether two objects (binary-single or binary-binary) will 
interact strongly. These strong interactions are performed 
using either simple recipes based on cross sections, or a 
dynamical integrator. For most of the work reported here, 
we use cross sections for the treatment of close binary- 
binary and binary-single interactions. These cross sec- 
tions were compiled from analytic fits to the results of 
numerical scattering experiments available in the litera- 
ture. Given the very large parameter space, reliable cross 
sections are available only for equal-mass encounters, and 
so we study only single-component clusters in this paper. 
All single stars arc assumed to have the same mass, and 
all binaries contain two identical stars with the same mass 
as the background single stars. All stars are treated as 
point masses, i.e., we neglect physical collisions between 
stars during interactions (cf. Bacon et al. 1996; Cheung 
et al. 2003). Our implementation follows closely that used 
in the Fokker-Planck study by GGCM91, which will serve 
as the main comparison for our work. 

2.1. Units and Definitions 

In our code we use the system of units defined by setting 
G = Mq = —4Eq = 1, where Mo is the initial cluster mass, 
and E is the initial cluster energy (excluding the binding 
energy in binaries). The corresponding unit of time is 

then idyn(0) = GM„ /2 (-4E )~ 3/2 • However, the natural 
timescale for cluster evolution is the relaxation time 

N 



When reporting results, however, we scale all times to the 
initial half-mass relaxation time, which we calculate using 
the standard definition given by Spitzer & Hart (1971), 



*r(0) 



-tdyn(0) . 



(1) 
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where 7 is a constant of order unity that must be deter- 
mined experimentally This relaxation time is used as the 
time unit in the Monte Carlo code. Therefore the ln(7^Vo) 
dependence factors out from all expressions used in simu- 
lating two-body relaxation (see Henon 1971b, or Paper I). 



*rh(0) 
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(2) 

where rh(0) is the initial half-mass radius of the clus- 
ter. Since t dyD (0) ~ [r£(0)/ (GM )} 1/2 , we have i rh (0) - 
0.1t r (0), where the numerical coefficient depends on the 
value of rh(0) for the particular initial model. 

When calculating rates for processes that do not occur 
on the relaxation timescale, such as dynamical interac- 
tions, one must adopt a specific value for 7 in the Coulomb 
logarithm when converting between dynamical and relax- 
ation times. In all our simulations for this paper we use 
7 = 0.4, the standard value adopted in most previous 
Fokker-Planck simulations, including those of GGCM91 
(cf. Paper I, where we show that 7 ~ 0.1 provides the best 
agreement with direct TV-body simulations for the evolu- 
tion of a single-component cluster to core collapse). In ad- 
dition, in calculating the binary-binary and binary-single 
interaction rates, GGCM91 use N c , the current number of 
stars in the cluster core, instead of No in the denominator 
of eq. 1 when converting between dynamical and relax- 
ation times. Although there is no rigorous justification for 
this choice, it appears reasonable, since the interactions 
occur mainly inside the high-density core, and we adopt 
the same prescription in our simulations. 

To estimate core quantities, including the number of 
stars in the core, we first sample over a small number 
of stars, typically 0.1-1% of the total number of stars in 
the cluster, to calculate the central density, po, and veloc- 
ity dispersion, cr . Since the velocity dispersion varies so 
slowly away from the center, we estimate the core velocity 
dispersion as cr c ~ er - The core radius is then defined to 
be 

and the number of stars in the core is calculated as 
4^ = 2^ 

3m 3m ' w 

where in is the average mass of a star in the cluster, and 
p c ~ 0.5p (Spitzer 1987). 

2.2. Binary-Single Interactions 

In a single time step, the probability that a binary will 
strongly interact with another object (single or binary) is 
given by 

P = awnAt , (5) 

where a is the cross section for the interaction, w is the 
relative velocity at infinity, n is the local number den- 
sity of stars (single or binary), and At is the time step. 
For binary-single interactions, a — <7b s , the binary-single 
interaction cross section, and n — n s , the local number 
density of single stars. In our code we calculate n s using 
a local sampling procedure, and take w to be the relative 
velocity between the nearest single star and the binary. 
The total cross section for close binary-single interactions 
is computed as ab s = ^Lx- Here 6 max is the impact pa- 
rameter which gives a distance of closest approach between 
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the binary and the single star of r m i n = 3.5 a, where a is 
the binary semi-major axis. For a binary of mass raj and 
single star of mass m we have 

2G(m + m b ) 



b z = r z ■ 1 + 

max min I 1 



w 2 r min 



(6) 



The coefficient of 3.5 is chosen such that all interactions 
with a distance of closest approach greater than r m i n result 
in only negligible energy transfer from the binary to the 
passing star in a fly by (see, e.g., Heggie 1975). As long as 
it is sufficiently large, the precise value of the coefficient 
has very little influence on the results. 

The binary-single interaction is performed only if a uni- 
form random number between and f is less than the 
computed probability. The interaction can in principle be 
computed exactly using a three-body dynamical integra- 
tor; this approach has many benefits, especially in provid- 
ing an accurate way of distinguishing between the various 
possible outcomes (see below). However, it also requires 
significantly more computational resources than using a 
simple analytic prescription. Following GGCM91, in these 
equal-mass simulations, we assume that the only outcome 
is binary hardening, and we use a semi-analytic fit (Spitzer 
1987, eq. [6-27]) to numerical results (Hut 1984) to com- 
pute the translational energy released. Let y = Ae/e be 
the fraction of the binding energy of the binary that is 
released as translational energy. The differential cross sec- 
tion for the interaction is given by (Gao et al. 1991), 



dy 



12.48™^ 



-) (i + y)-V - 5 . 



(7) 



where the critical velocity is w cr = (3Gm/2a) 1 / 2 , the ve- 
locity at infinity for which the total energy of the system 
is zero and complete ionization is possible. The quan- 
tity y is drawn randomly from eq. (7) using the rejection 
technique. With this recoil energy and a scattering angle 
drawn at random in the center-of-mass frame, new veloc- 
ities and orbital energies in the cluster are calculated for 
the emerging binary and single star. 

As a first step toward the eventual goal of treat- 
ing all binary interactions exactly, we have incorporated 
into our code a dynamical integrator to perform binary- 
single interactions. Specifically, we use the three-body 
integrator scatter3 from the Starlab software environ- 
ment (see Appendix B of Portegies Zwart et al. 2001, 
and http://www.manybody.org). Scatter3 uses a time- 
symmetrized Hermite integrator and analytical continua- 
tion of unperturbed orbits to evolve the three-body system 
until an unambiguous outcome is obtained. The main ben- 
efit of using an exact treatment is the ability to study non- 
equal-mass systems, although for comparison with cross 
sections we restrict ourselves to the equal-mass case here. 
The implementation of the three-body integrations follows 
that of cross sections: first, the probability for an en- 
counter to occur is calculated according to eqs. (5) and 
(6); next, with velocity at infinity w, the impact param- 
eter b is chosen randomly in area, i.e., with probability 
dP{b) = 2irb db / (7r6^ ax ) . The binary eccentricity is as- 
sumed to follow a thermal distribution with dP(e) = 2e, 
and all angles are chosen assuming random orientation and 
phase. The dynamical interaction is then calculated and 
its outcome is used to determine the new binding energy 
of the binary and the new orbits for the binary and single 



star in the cluster. The current implementation properly 
handles the outcomes of preservation and exchange, but, 
for the sake of comparison with cross sections, currently 
ignores ionizations. This is justified here since we con- 
sider only hard primordial binaries in our simulations (see 
Sec. 3.1)'. 

2.3. Binary-Binary Interactions 

To calculate the probability that a close binary-binary 
interaction should occur in a time step, we use eq. (5) with 
c = 0bb, the binary-binary interaction cross section, and 
n = rib, the local number density of binaries. We calcu- 
late nb using a local sampling procedure, and we take w to 
be the relative velocity between the current binary and the 
nearest binary. Following GGCM91, for the binary-binary 
interaction cross section we use the results of Mikkola 
(1983a, b, 1984a, b) for encounters between equal-mass bi- 
naries. In the case where one binary has a much higher 
binding energy than the other (ei ;§> £2), Mikkola (1984a) 
provides a semi-analytic fit to his numerical results, giving 
a collision cross section 



Obb w 16.6 



In 



29|e 2 



mw 2 + 0.04|e 2 | 



2/3 



Gma--, 



(8) 



where w is the relative velocity of the two binaries at in- 
finity, m is the mass of each star in the binaries, and 02 is 
the semi-major axis of the softer binary. 

An interaction between two binaries can result in many 
possible outcomes. Since we consider only hard binaries 
in this study, the most probable outcomes are (1) disrup- 
tion of the softer binary and hardening of the harder one, 
and (2) the formation of a stable hierarchical triple with 
a single star ejected. As much as ~ 1/3 of close binary- 
binary encounters may result in the formation of such a 
triple. However, for a triple system to remain long-lived 
in the dense cluster environment, its outer orbit must be 
sufficiently tight. The formation of a long-lived hierar- 
chical triple is expected to be much less common in a 
dense cluster (see Ford et al. 2000). Therefore we assume 
for simplicity that all hierarchical triples formed are dis- 
rupted immediately. The only outcome of binary-binary 
encounters that we treat in our simulations is then case 
(1) above. Mikkola finds that, on average, approximately 
one half of the combined binding energy (ei + £2) of the 
two binaries is released in the form of translational energy, 
AE t = y(e\ + e 2 ). The semi-analytic fit given by Mikkola 
(1984a) for the distribution of translational energies pro- 
duced is rather complicated. Instead, we use a simplified 
version of the distribution, 



G(y) 



49 



-y 



-11/4 



(9) 



proposed by GGCM91. The mean value of this distribu- 
tion, (y) ~ 0.47, is in good agreement with the results 
of Mikkola for interactions resulting in a binary and two 
single stars. 

We also adopt a simplified overall binary-binary colli- 
sion cross section, by replacing the expression in square 
brackets in eq. (8) by its value at e 2 = jmw 2 , yielding 



ebb = 31.8- 



Gmao 



(10) 
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The energy required to disrupt the softer binary, as well as 
the total translational energy released in the collision AE t , 
are both generated at the expense of the surviving binary. 
Thus the binding energy of the surviving pair increases by 
an amount 62 + v{^i + £2)- According to Mikkola (1983a), 
for collisions between binaries of equal binding energies 
producing a binary and two single stars, typically about 
1/4 of the translational energy produced is carried away 
by the binary, and the remaining is distributed randomly 
among the two single stars. For simplicity, we assume that 
this prescription is applicable to collisions between bina- 
ries of unequal binding energies as well. We select the 
direction of the recoil velocity between the binary and the 
single stars randomly in the center-of-mass frame. 

If a binary does not undergo a strong interaction with a 
single star or another binary, it is then treated as a single 
star in the usual two-body relaxation step (see Paper I), 
during which its internal structure is left unchanged. 

3. RESULTS 

3.1. Initial Conditions and Summary of Model Results 

For our initial cluster models we use both the Plummcr 
model, assumed to be isolated (i.e., with no tidal boundary 
enforced), and a variety of tidally truncated King models, 
assumed to be on a circular orbit in the Galaxy (i.e., with 
a fixed external tidal potential). Mass loss through the 
tidal boundary is treated as in Paper II, using a criterion 
based on the apocenter distance of each stellar orbit in 
the cluster, and an iterative procedure to determine both 
the mass loss and the new position of the tidal boundary 
after each relaxation timestep. The initial binary fraction 
/b (defined as the fraction of stars, by number, that are 
binaries) varies between and 30%. In a few cases, for 
calibration, we have also performed simulations in which 
the binaries are present, but all interactions are turned 
off; these models are equivalent to two-component mod- 
els in which a small fraction of (single) stars have twice 
the mass of the background stars (see Watters et al. 2000 
and Fregeau et al. 2002 for other studies of two-component 
clusters using our code). 

The binaries are distributed initially in the cluster ac- 
cording to the same density profile as for single stars. 
Hence no initial mass segregation is assumed for the bina- 
ries. The distribution of the internal binding energy of the 
binaries is assumed to be uniform in log e between a min- 
imum value e m i n and a maximum value e max . Following 
GGCM91, we consider only hard binaries, with the min- 
imum binding energy e m i n = mer c (0) 2 , where ct c (0) is the 
initial central velocity dispersion. Soft binaries, if present, 
would be assumed to be ionized (destroyed) as soon as they 
participate in a strong interaction. Therefore they would 
not affect the overall evolution of the cluster significantly. 
For the maximum binding energy we take 6 max — 133e m i n , 
which is approximately the binding energy of a contact 
binary for two solar-like stars if <t c (0) ~ 10 kms -1 . The 
precise value of e max has little influence on our results, 
since very hard binaries behave essentially as single more 
massive stars (with very small interaction cross section). 

Table 1 lists the parameters of the main models we con- 
sidered, as well as the main results of our simulations for 
each cluster. The first column identifies the initial clus- 
ter model, Plummer or King, and the value of the con- 



centration parameter Wo (dimensionless central potential) 
for King models. The second column gives the initial bi- 
nary fraction /b. All simulations were performed with 
N = 3 x 10 5 stars (including binaries) initially in the clus- 
ter. The following columns summarize the main results of 
our dynamical simulations. For each model we first give 
the time of core collapse t cc , in units of the initial half- 
mass relaxation time i r h(0), defined by cq. (2). Here core 
collapse is defined as the moment when the core density 
reaches its first maximum. This can be determined typi- 
cally to within a statistical error of at most a few percent 
in our simulations. We then give the total mass of the 
cluster at the moment of core collapse (in units of its ini- 
tial total mass), and the fraction of binaries that remain 
at that moment. For clusters that disrupt completely be- 
fore reaching core collapse, we list the disruption time idis 
instead of t cc . 

3.2. Comparison with Direct Three-Body Integration 

As a simple test of our code and the approximate treat- 
ment of interactions, we compare the use of cross sections 
with dynamical integrations of binary-single encounters. 
In future work, we will also implement dynamical integra- 
tions of binary-binary interactions, and we will use more 
detailed comparisons to re-calibrate the various recipes 
based on cross sections. Here our intent is merely to 
demonstrate that these simple recipes are reasonably accu- 
rate. We have not changed our prescriptions to try to bet- 
ter match the results of the dynamical integrations, since 
a main goal in this first study is to provide comparisons 
with the Fokker-Planck simulations of GGCM91 that used 
the same simple prescriptions. For this test binary-binary 
interactions were turned off. In reality, they tend to dom- 
inate the energy production (see Sec. 3.3). Thus this sim- 
ple test also allows us to study specifically the effects of 
three-body interactions on the overall cluster evolution. 

Figure 1 shows the evolution of an isolated cluster de- 
scribed initially by a Plummer model with N = 3 x 10 5 
stars and 20% binaries. Solid lines correspond to the sim- 
ulation using direct three-body integrations, while dashed 
lines show the results using our simple cross sections. The 
top panel shows the total mass in binaries in the clus- 
ter, decreasing as binary burning proceeds. Since all bi- 
naries in the model are hard, binary-single interactions 
(unlike binary-binary interactions) cannot destroy a bi- 
nary, and therefore binaries can only be lost by ejection 
from from the cluster (typically following significant hard- 
ening through multiple interactions; see Hut et al. (1992b) 
and Sec. 3.5 below). The rate of binary ejection acceler- 
ates abruptly at i/i r h(0) ~ 8 — 10 near core collapse. The 
middle panel of Figure 1 shows the energy generated in 
binary-single interactions, as a fraction of the total initial 
binding energy of the cluster. By the time of core collapse, 
this is only ^0.1. This amount of energy is not sufficient 
to delay core collapse significantly. In fact the binaries, 
through mass segregation, accelerate core collapse in this 
(artificial) simulation (recall that the core collapse time 
of a single-component Plummcr model without binaries is 
given by i cc Arh(0) — 14). The bottom panel shows vari- 
ous characteristic radii in the cluster: from top to bottom, 
the half-mass radius of single stars, the half-mass radius 
of binaries, and the core radius, all in units of the initial 
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half-mass radius. 

The agreement between the two methods is strong, al- 
though the total energy generated in binary-single interac- 
tions is slightly smaller when calculated by direct dynami- 
cal integrations. Consequently, the model using dynamical 
integrations reaches core collapse sooner than the model 
using cross sections, because less energy is generated to 
support the core against collapse. We believe that this 
difference comes from the deterministic treatment of bi- 
nary hardening with cross sections, in which every binary- 
single interaction results in a hardened binary. In reality 
the widest hard binaries in the simulation, which are right 
around the hard/soft boundary (and have the largest in- 
teraction cross section) have roughly equal probabilities of 
hardening and softening in an interaction (Heggie 1975). 
To partly restore consistency between the two treatments, 
we ignore dynamical integration outcomes that result in 
ionization of the binary. Were these included, the total 
energy generated in binary-single interactions would de- 
crease further, by roughly 50%. This would cause the bi- 
nary population to become more centrally concentrated. 
Thus, in a realistic cluster simulation, we would expect 
the ratio of the energy generated by binary-single inter- 
actions to binary-binary interactions to decrease by more 
than 50% compared to predictions of cross-section based 
recipes. 

3.3. Isolated Clusters 

We consider first the evolution of Plummcr models con- 
taining N = 3x 10 5 stars with a range of binary fractions 
/b- As a further test of our method, we show in Fig. 2 the 
evolution of the various energies and the virial ratio of a 
system with /b = 0.1. Since dynamical relaxation is not 
built into our numerical method, the degree to which virial 
equilibrium is maintained during a simulation is our most 
important indicator of numerical accuracy. We monitor 
this, as well as energy conservation, in all our runs, and 
terminate a calculation whenever these quantities deviate 
from their expected values by more than a few percent 
(this typically happens when the number of binaries has 
been reduced to a very small value, or, in tidally truncated 
clusters, when the total number of stars remaining in the 
cluster becomes very small; See Sec. 3.4 below). 

Figures 3, 4 and 5 show the evolution of models with 
fb = 0.02, 0.1, and 0.2, respectively. The main impact of 
introducing binaries in the models is very clear: core col- 
lapse is delayed considerably. Even for a cluster with only 
2% binaries initially (Fig. 3), t cc increases by more than a 
factor 2. Clusters with Jb — 0.1 — 0.2 can avoid core col- 
lapse for ~ 100 t r h(0) (Figs. 4 and 5). For the vast majority 
of globular clusters in our Galaxy, where i r h(0) ~ 10 9 yr, 
this timescale exceeds a Hubble time. If all globular clus- 
ters in our Galaxy were born with /b > 0.1, only those 
with very short initial relaxation times would have had a 
chance to reach core collapse. However, for real clusters, 
tidal truncation and mass loss (Sec. 3.4) as well as stellar 
evolution (Paper II) complicate this picture considerably. 

In Figure 3, we also show for comparison the evolu- 
tion of the core radius for a model in which binaries are 
present but all interactions are turned off (short-dashed 
line in the bottom panel). Even with a binary fraction 
as small as /b = 0.02 in this case, core collapse occurs 



significantly earlier than in a single component Plummcr 
model (at £ C cArh(0) — 10 instead of 14). This shows the 
expected tendency for the heavier component of binaries 
to accelerate the evolution to core collapse, and the result 
is in good agreement with previous studies of core col- 
lapse in two-component clusters (see, e.g., Watters ct al. 
2000). Note that for sufficiently large binary fractions, 
these two-component models become "Spitzer-unstable," 
i.e., the core collapse is driven entirely by the heavier com- 
ponent. Using the stability criterion derived by Watters 
et al., A = (M 2 /M 1 )(m 2 /m 1 ) 2A < 0.32, here with an indi- 
vidual mass ratio m^/mi — 2 and a total component mass 
ratio M2/M1 = 2/b/(l — fb) we expect the Spitzer insta- 
bility to appear whenever the binary fraction /b > 0.03. 
Thus all our models with binary fractions above a few 
percent should evolve on a relaxation timescale to a state 
where the dynamics of the cluster core is largely dominated 
by the binaries. Indeed, looking at the middle panels of 
Figures 4 and 5, we see that, with /b = 0.1 — 0.2, the 
energy generation is largely dominated by binary-binary 
interactions. In contrast, for /b = 0.02 (middle panel of 
Fig. 3), binary-binary and binary-single interactions con- 
tribute roughly equally. 

To quantify the effect of primordial binary burning on 
the core collapse time, we have repeated calculations with 
binaries present but all interactions turned off for six differ- 
ent models with varying /b. For each model, we can then 
properly calculate the ratio of core collapse times with 
and without binary interactions (but with mass segrega- 
tion effects present in both cases). The results are shown 
in Figure 6, where this ratio is plotted as a function of the 
binary fraction. A simple linear fit gives 

t cc ~t cc (f b = 0) x (75/ b + l), 

and reproduces the numerical results to within ~ 30% 
in the range /b = — 0.3. The notation we use here, 
l %c{fb = 0)," means the core collapse time of a cluster 
with the same fraction of "inactive" binaries, rather than 
with no binaries. 

Also shown in Figure 3 for comparison is the result of 
a simulation in which all interactions are included, but 
binary-single interactions are calculated by direct three- 
body integrations as in Sec. 3.2 (long-dashed line in the 
bottom panel). This comparison is useful again as a test of 
the simple treatment based on cross sections, since binary- 
single interactions play an important role as a source of 
energy in this model with /b = 0.02. Our conclusion is 
the same as in Sec. 3.2: the agreement is very good until 
i/i r h(0) ~ 15, but then the two simulations diverge and 
the model computed with direct three-body integrations 
collapses slightly earlier (t cc /irh(0) — 19 instead of 22). 
In spite of this slight offset, after core collapse the core 
re-expansion and gravothermal oscillations also look very 
similar in the two simulations. 

The top panels in Figures 3-5 show the evolution of 
the total cluster mass, as a fraction of the initial mass, 
and the remaining mass in binaries, as a fraction of the 
initial mass in binaries (this is also the remaining frac- 
tion by number since all binaries have the same mass). 
Binary-binary interactions are the main process respon- 
sible for the destruction of binaries in these simulations 
(since the softer binary is assumed to be disrupted in each 
interaction) . In the absence of evaporation through a tidal 



MONTE CARLO SIMULATIONS OF GLOBULAR CLUSTER EVOLUTION. III. 



7 



boundary, mass loss from the cluster comes almost entirely 
from stars and binaries ejected through recoil following an 
interaction. The mass loss rate therefore increases with in- 
creasing binary fraction. At core collapse, the total mass 
loss fraction is about 5%, 15%, and 25% for /b = 0.02, 
0.1, and 0.2, respectively. However, while the total num- 
ber of binaries in the cluster continuously decreases, the 
remaining fraction of binaries at core collapse appears to 
be roughly constant, around 0.2, independent of /b- This 
is sufficient to power many cycles of gravothermal oscil- 
lations after the initial core collapse, even for initial bi- 
nary fractions as small as a few percent. For ,/b = 0.2 
(Fig. 5), we were able to extend our numerical integration 
all the way to almost ~ 10 3 i r h(0), at which point several 
thousand binaries are still present in the central region of 
the cluster (this timescale would of course vastly exceed a 
Hubble time for most Galactic globular clusters!). 

The bottom panels in Figures 3-5 show the evolution of 
several characteristic radii. The most important is the core 
radius r c (recall that, by our definition, cq. (3), the central 
density scales approximately as po oc r c ~ 2 since the cen- 
tral velocity dispersion is approximately constant). Even 
in deep core collapse, the core radius of our models never 
decreases by more than a factor ^100 (corresponding to 
an increase in the central density by ~ 10 4 ). Models with 
higher binary fractions contract very little (see Fig. 5: the 
first and deepest core collapse corresponds to a decrease in 
r c by less than a factor 10). Also shown are the half- mass 
radii of the binaries 7"h,b and of the single stars r"h, s - The 
half-mass radius of the single stars always increases mono- 
tonically for these isolated clusters. In contrast, the half- 
mass radius of the binaries tends to increase on average 
but shows a much more complex behavior that depends 
strongly on the binary fraction and on the particular dy- 
namical phase in the evolution of the cluster. The trend 
is for rh,b to decrease during normal cluster evolution, as 
the binaries mass segregate to the cluster core, and to in- 
crease dramatically during core collapse, as the density of 
binaries in the core grows and the rate of binary-binary 
interactions grows more quickly than the rate of binary- 
single interactions. This causes many softer binaries in 
the core to be disrupted and many harder binaries to be 
ejected out of the core through recoil. For sufficiently low 
binary fractions (Figs. 3 and 4), rh,b eventually becomes 
larger than rh, s after core collapse. For high binary frac- 
tions (Fig. 5), the binaries remain always much closer to 
the center of the cluster. 

We now turn to a more detailed discussion of the Plum- 
mer model with 10% binaries (Fig. 4), including a com- 
parison with the Fokker-Planck results of GGCM91 (see 
their Figs. 1-3), who consider this their "standard model." 
Qualitatively, our results are in very good agreement up to 
core collapse. After an initial phase of contraction lasting 
~ 10t r h, the core radius becomes nearly constant and the 
cluster enters a long phase of quasi-thermal-equilibrium. 
This is the stable "binary burning" phase, analogous to 
the main sequence for a star. During this phase, the rate 
of energy production through interactions in the core is 
balanced by the rate at which energy flows out in the 
outer halo, which continuously expands (in the absence of 
a tidal boundary). Core collapse occurs rather suddenly 
at the end of this phase. GGCM91 find i cc /t rh (0) ~ 50 



for this model, while we find t cc /t T h(0) — 70. This ini- 
tial core collapse is followed by gravothermal oscillations, 
which are clearly still powered by primordial binary burn- 
ing. We were able to follow these oscillations accurately 
until t/*rh(0) > 200, while GGCM91 terminate their cal- 
culation at i/i r h(0) — 90. 

Upon closer examination and quantitative comparisons, 
some more significant differences become apparent. First, 
we see that the initial contraction phase appears much 
deeper in the model of GGCM91, with r c decreasing by 
almost an order of magnitude, while in our model the core 
contracts by a factor of about 3. During the stable binary 
burning phase, the cluster also appears somewhat more 
centrally concentrated in the model of GGCM91. Just 
before core collapse, they find r c /?*h(s) ~ 0.01, while our 
value is ~ 0.04. On the other hand, the rate of binary 
burning and destruction is nearly the same in the two mod- 
els. Compare, for example, the evolution of Mb/Mb(0) 
in Figure 4 to the same quantity plotted in Fig. 2a of 
GGCM91. Although there are slight differences in the 
shapes of the two curves, the reduction to 0.8 occurs after 
about 10 i r h(0) in both cases, and the reduction to 0.5 after 
about 28t r h(0). By i/i r h(0) ~ 70 the number of binaries 
has been reduced to 0.2 of its initial value in both models. 
This agreement is especially surprising since in our model 
this is still (just) before core collapse, while in GGCM91's 
model several cycles of gravothermal oscillations have al- 
ready occurred. 

There are several reasons to expect differences between 
our results and those of GGCM91's Fokker-Planck simu- 
lations, even though our treatments of individual binary- 
single and binary-binary interactions are essentially iden- 
tical. 

First, GGCM91's representation of binaries is in terms 
of a separable continuous distribution function in E, the 
orbital energy in the cluster, and e, the internal energy of 
the binary. In fact, there is a strong and complex correla- 
tion between a binary's binding energy and its position in 
the cluster (or cquivalently its energy in the cluster), with 
harder binaries concentrated near the cluster core (see Hut 
et al. 1992b and Sec. 3.5). We suspect that GGCM91's 
choice of a separable distribution function has the effect 
of reducing the energy generation rate, since then pro- 
portionately more soft binaries will be chosen for binary 
interactions — interactions that predominantly liberate a 
constant fraction of the total binding energy available (see 
Sec. 2 and Heggie 1975). 

Second, 1-D Fokker-Planck results are known to differ 
from 2-D results in general (most notably in the predic- 
tion of the rate of tidal stripping; see Paper II) . Enforcing 
isotropy in the stellar velocity distribution is likely to af- 
fect the dynamics of the core around the time of collapse, 
when this distribution may be changing rapidly and the 
increased interaction rate may be causing anisotropy. In- 
deed, Baumgardt et al. (2003) have recently found with 
A^-body simulations that anisotropy near the cluster cen- 
ter becomes significant during core collapse. 

Third, the only explicit dependence on N in the Fokker- 
Planck approach is through the Coulomb logarithm, so, 
even though GGCM91 set N = 3 x 10 5 for their treatment 
of interactions, it is not clear in what sense their results, 
which assume a smooth, continuous distribution function, 
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correspond to this particular value of A. 

Finally, we point out that our results for the initial con- 
traction and core size during the binary-burning phase are 
in much better agreement with those of direct A-body sim- 
ulations including primordial binaries. Indeed, both Hcg- 
gie & Aarseth (1992) (see their Figs. 5 and 18) and McMil- 
lan et al. (1990) (see their Fig. 1) find, as we do, that the 
cluster core contracts typically by a factor of about three 
from its initial size. GGCM91, on the other hand, find core 
contraction by an order of magnitude, a direct indication 
that their method underestimates the energy generation 
rate. 

Perhaps a more significant difference between our re- 
sults and those of GGCM91 is in the post-collapse evo- 
lution. GGCM91 find much more frequent, erratic, and 
deeper gravothermal oscillation cycles. Our model shows 
almost quasi-periodic oscillations with period ~ 40 i r h and 
peak-to-peak amplitude r C!max /r C!m j n ~ 100. Instead, 
GGCM91 find 7 oscillations of widely varying periods be- 
tween t/t r h(0) = 50 and 90, and r c , max /r Cjm i n ~ 10 3 . We 
believe this may possibly be due to differences in the nu- 
merical method of calculating N c (for use in eq. (1) — see 
discussion in Sec. 2.1), although it is not clear from their 
paper what method they actually use. To check the va- 
lidity of our results near core collapse, we have examined 
more carefully the dynamics governing core re-expansion 
after collapse. In Figure 7, we show the evolution of the 
temperature profile in the cluster as the system undergoes 
a core collapse and rebound. The temperature in the clus- 
ter normally decreases outward everywhere, as in a star in 
thermal equilibrium. However, during deep core collapse, 
a "temperature inversion" develops for a short time. This 
temperature inversion is responsible for driving the rapid 
re-expansion of the core, as energy is now flowing inward. 
This mechanism has been predicted theoretically for a long 
time (Sugimoto & Bettwieser 1983; Heggie & Ramamani 
1989), and has been observed directly in recent A-body 
simulations (Makino 1996). However, to our knowledge, 
ours is the first numerical demonstration of this effect for 
a cluster containing primordial binaries. In all previous 
studies, the binaries were assumed to form dynamically via 
three-body interactions or two-body tidal captures during 
deep core collapse. As noted in the introduction, these 
mechanism are now considered unrealistic, as they most 
likely lead to stellar mergers (which were not taken into 
account in the previous studies). 

3.4. Tidally Truncated Clusters 

We now present our results for more realistic, tidally 
truncated clusters. Figures 8-11 show the evolution of 
King models with Wo — 7 and initial binary fractions 
/b = 0.02, 0.05, 0.1, and 0.2, respectively. Several striking 
differences with isolated clusters are immediately appar- 
ent. First, we see that the initial core contraction phase 
is absent. For /b < 0.1, the core radius decreases slowly 
and monotonically all the way to collapse. For higher bi- 
nary fractions (Fig. 11), the core radius increases slightly 
at first. This is simply because the initial binary burning 
rate in this model is close to the rate needed to reach ther- 



mal equilibrium. The higher the initial binary fraction and 
central density (see below), the stronger the tendency for 
the core to expand initially instead of contracting. Second, 
core collapse 5 or complete disruption always occurs in less 
than about 45t rn (0). For /b > 0.1, the disruption time tdis 
decreases with increasing binary fraction, and complete 
disruption occurs before any deep core collapse. This is in 
contrast to the models with lower binary fractions (Figs. 8 
and 9), where core collapse followed by gravothermal os- 
cillations (similar to those observed for isolated clusters in 
the previous section) occur before disruption. The pos- 
sibility for a cluster to suffer complete disruption before 
core collapse is a qualitatively new behavior introduced 
by primordial binaries. Indeed, all King models without 
binaries (and without stellar evolution) reach core collapse 
before disrupting (see Paper II and Quinlan 1996). 

Figures 12 and 13 show the evolution of King models 
with 10% binaries but different values of Wq. For the very 
centrally concentrated cluster with Wo = 11 (Fig. 12), sig- 
nificant core expansion occurs in the first few relaxation 
times (with r c increasing by about an order of magnitude) . 
This is a more extreme example of the behavior already 
noted in Fig. 11. The final evolution of this cluster is 
also peculiar: this is one of few examples (see Table 1) we 
encountered where the binaries are completely exhausted 
before the cluster disrupts. At t/t T h{0) ~ 30, about 20% 
of the initial cluster mass remains in single stars, and the 
cluster undergoes deep core collapse. Since there are no bi- 
naries left, and our simulations include no other source of 
energy, no re-expansion can occur and we must terminate 
the calculation. 

For a model with Wo = 3, which has a much more 
nearly-uniform density profile initially, complete disrup- 
tion occurs before core collapse even with much lower bi- 
nary fractions (see Table 1). For the model with /b = 0.1 
shown in Fig. 13, disruption occurs at i/i r h(0) — 15. For 
comparison a single-component Wo = 3 King model un- 
dergoes core collapse at i/i rn (0) ~ 12 (Paper II). Also 
note how the core contracts throughout the evolution at 
a nearly-constant rate much faster than in models with 
higher values of Wo (compare, e.g., Fig. 10). Just before 
final disruption, the core radius has decreased by about an 
order of magnitude from its initial value. 

Only a few small A-body simulations of tidally trun- 
cated clusters with primordial binaries have been reported 
previously (Heggie & Aarseth 1992; McMillan & Hut 
1994). Detailed comparisons are not possible because 
these studies assumed rather different initial models and 
the A-body results (for A ~ 1000 — 2000) are very noisy. 
However, we do see good qualitative agreement, with mass 
loss rates ~ 10 times larger than for isolated clusters, and 
complete disruption also observed after a few tens of ini- 
tial half-mass relaxation times in all A-body simulations. 
We are not aware of any previous Fokker-Planck simula- 
tions of tidally truncated clusters with primordial binaries 
(GGCM91 considered only isolated Plummer models). 

3.5. Evolution of the Binary Population 



5 It should be noted that, in much of the literature on globular cluster dynamics, the term "core collapse" is used to refer to the initial 
core contraction phase (which we see here can actually be expansion instead), and what we call the binary burning phase is then called the 
"post-collapse" phase. Clearly this terminology no longer makes sense, and should be abandoned. What we call "core collapse" in this paper 
refers to the brief episodes of deep core collapse at the onset of and during gravothermal oscillations. 
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In addition to affecting the global cluster evolution, bi- 
nary interactions also strongly affect the properties of the 
binaries themselves. The study of the evolution of a pri- 
mordial binary population dates back to the seminal work 
of Heggie (1975), but it is only recently that detailed nu- 
merical simulations of large binary populations in globular 
clusters have been performed (Hut et al. 1992b; Giersz & 
Spurzcm 2000). We can use the results of our Monte Carlo 
simulations to study the dynamical evolution of the binary 
population. 

Fig. 14 shows the evolution of the binary fraction, /b, 
in different regions of our evolving King models. There is 
a clear trend for /b to increase in the core and decrease 
in the halo with time, as well as a trend for /b to grow 
more with smaller Wq- In spite of mass segregation and 
the tendency for binaries to dominate the central region 
of a cluster (following the development of the Spitzer in- 
stability; see Sec. 3.3), core binary fractions rarely exceed 
0.5 in our models. Thus the range of initial binary frac- 
tions we consider are at least in rough agreement with the 
measurements of core binary fractions in globular clusters 
today (~ 0.1 — 0.4; see Sec. 1). Note, however, that the 
present-day binary fraction in the core of a cluster cannot 
be related simply to the cluster's initial binary fraction, 
as it may depend in a complicated way on several initial 
parameters. For example, we see that an initial Wo = 11 
model with /b = 0.1 has, during most of its evolution, 
about the same core binary fraction (~ 0.15) as an ini- 
tial Wq = 7 model with /b = 0.05. In addition, recall 
that our definition of /b in these simulations includes only 
the hard binaries. For reasonable distributions of primor- 
dial binary separations, including several more decades on 
the soft side, the true initial binary fraction in the cluster 
might have been ^2 — 3 times larger than our quoted value 
of /b (this is the reason why we did not consider values 
of /b > 0.3, which could not be realistic, unless dynamics 
already plays an important role during the process of star 
formation; see Clarke et al. 2000). 

Fig. 15 shows the evolution of the primordial binary 
population in a Wo = 7 King model with /b = 0.2. Each 
2-D histogram shows the distribution of binding energies 
(initially flat in log e) and radial positions in the cluster. 
In addition to the clear tendency for mass segregation and 
hardening of the binaries, we note the development of a 
strong correlation between hardness and radial distribu- 
tion: harder binaries tend to concentrate in the cluster core 
much more than softer binaries, in spite of having all the 
same mass (and in contrast to the fundamental assumption 
made in the Fokker-Planck calculations of GGCM91). The 
general trends observed here are in good qualitative agree- 
ment with the results of previous studies using more ide- 
alized models (Hut et al. 1992b; Giersz & Spurzem 2000, 
see, e.g., their Fig. 25). 

Near the end of the evolution shown in Fig. 15 (but 
already ~ 10i r h(0) before complete disruption, when the 
cluster still retains about 40% of its initial mass), a partic- 
ularly striking situation develops where all the surviving 
binaries in the cluster core are extremely hard. Recall that 
our initial upper limit on the binding energy of a binary 
(~ 10 2 kT, where "fcT" is the average kinetic energy of 
stars in the core) roughly corresponds to contact for two 
solar-like stars. Therefore, most of the binaries remaining 



after ~ 30t r h(0), with binding energies now in the range 
~ 10 2 — 10 3 kT, would have merged if they contained solar- 
like stars (perhaps forming blue stragglers). Of course, in 
a real cluster, many of these binaries could contain com- 
pact objects (most likely heavy white dwarfs and neutron 
stars) and would then have survived. We cannot address 
any of these issues here, since our simulations are clearly 
too idealized, but we point out that globular cluster cores 
are indeed observed to contain large populations of blue 
stragglers, WUMa binaries (eclipsing systems containing 
two main-sequence stars in a contact configuration; see, 
e.g., Albrow et al. 2001), and a variety of "ultracompact" 
binaries containing neutron stars and white dwarfs (the 
most extreme example being perhaps the "11-minute" X- 
ray binary in NGC 6624; see, e.g., Deutsch et al. 2000). 

4. SUMMARY AND COMPARISON WITH OBSERVATIONS 

We have performed, for the first time, discrete simula- 
tions of globular clusters with realistic numbers of stars 
and primordial binaries, using our 2D Monte Carlo code 
with approximate analytical cross sections for primordial 
binary interactions. 

We have compared the use of cross sections with ex- 
act, dynamical integrations of binary-single encounters, 
and find that the agreement between the two methods is 
strong, although our current implementation of the cross 
sections, based on the Fokker-Planck study by Gao et al. 
(1991), tends to overestimate slightly the energy genera- 
tion rate. Consequently, models that use cross sections 
tend to overestimate core collapse times for clusters in 
which binary-single interactions dominate. However, we 
find that binary-binary interactions dominate the energy 
generation for /b > 0.03, a result that is in quantitative 
agreement with a simple Spitzer-type stability criterion 
applied to the component of binaries. 

We have studied the evolution of isolated clusters with 
varying binary fractions, and have found that the pres- 
ence of even a small fraction of binaries is sufficient to 
delay significantly the onset of core collapse. Isolated 
clusters with a binary fraction greater than about 0.1-0.2 
can avoid core collapse for as much as ~ 10 2 — 10 3 i r h(0). 
We find a simple linear relation between the core collapse 
time of a cluster, i cc , and the core collapse time of the 
same cluster with binary interactions turned off (but mass 
segregation still present), denoted t cc (/b = 0), given by 
t cc ~ i cc (./b = 0) x (75 /b + 1). We have compared our re- 
sults with those of Gao et al. (1991), and find reasonable 
agreement, with nearly identical rates of binary burning 
and destruction. Gao et al. (1991), however, find a shorter 
core collapse time, a deeper initial core contraction, and 
significantly more erratic behavior during the gravother- 
mal oscillation phase. We attribute the differences primar- 
ily to their neglect of the strong correlation between binary 
hardness and spatial distribution in the cluster, as well as 
fundamental differences between their 1-D Fokker-Planck 
method and our 2-D Monte Carlo method. Our results for 
the initial core contraction are in much better agreement 
with those of direct iV-body simulations. In addition, we 
have presented the first numerical demonstration of the 
theoretically predicted temperature inversion powering re- 
expansion after core collapse and gravothermal oscillations 
for a cluster with primordial binaries. 
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We have also considered more realistic, tidally truncated 
King models. We have found that the initial core con- 
traction phase is absent in these systems, or replaced by 
an initial expansion of the core, and that core collapse 
or complete tidal disruption always occurs in less than 
~ 50i r h(0). For a binary fraction > 0.1, the disruption 
time decreases with increasing binary fraction, and com- 
plete disruption occurs before any deep core collapse. The 
possibility for a cluster to suffer complete disruption before 
core collapse is a qualitatively new behavior introduced 
by primordial binaries. Our results are in good qualita- 
tive agreement with previous studies of tidally truncated 
clusters containing primordial binaries. 

We have already argued in Section 3.5 that our re- 
sults are in general agreement with current determina- 
tions of binary fractions in globular cluster cores, typically 
~ 0.1 — 0.4. We now briefly consider our basic predictions 
for the structural parameters of clusters during the binary 
burning phase and compare them to the observed struc- 
tural parameters of globular clusters. While our models 
are clearly far too idealized for any detailed comparison, 
it is useful to examine at least the most fundamental struc- 
tural parameters: the core radius r c , the half- mass radius 
rh, and, for tidally truncated clusters, the tidal radius r t . 
Since the overall scale is largely irrelevant (although it 
could in principle be set by relating the maximum bind- 
ing energy of the binaries to a stellar radius), we consider 
only the two ratios rh/r c and r t /r c (or, equivalently, the 
concentration parameter c = log(r t /r c ) often derived by 
observers using King model fits to photometric data). 

Fig. 16 shows the evolution of rh/r c and the concentra- 
tion parameter c for several King models (see Sec. 3.4). 
Fig. 17 shows distributions of rh/r c and c for Galactic 
globular clusters, with data taken from the compilation of 
Harris (1996). The top panel shows a histogram of rh/r c 
values for all Galactic globular clusters, including those 
classified observationally as "core-collapsed 6 ." The middle 
and bottom panels show the distributions of observed val- 
ues for rh/r c and c, respectively, with the "core-collapsed" 
clusters excluded. First, in Fig. 16, note the tendency for 
the concentration parameter in clusters with reasonable 



binary fractions (/b > 0.05) to remain around, or even 
converge to, c ~ 1.5. This is in reasonable agreement with 
the bottom panel of Fig. 17, which shows the observed dis- 
tribution also centered around c ~ 1.5. The top panel in 
Fig. 16 shows most clusters in the binary burning phase 
with rh/r c < 10, also in quite reasonable agreement with 
the observed distribution (Fig. 17 b), although the ob- 
served peak around rh/r c ~ 2 would require that most ini- 
tial models be less centrally concentrated than our Wo = 7 
King models. We also note that both c and rh/r c increase 
significantly, and sometimes dramatically, for clusters ap- 
proaching tidal disruption. Thus the suggestion from our 
results might be that clusters classified observationally as 
"core-collapsed" are those in the last few relaxation times 
before destruction in the Galactic tidal field. Most of our 
King models appear to spend roughly the last 10-20% of 
their lives with rh/r c > 10, or c > 2, again not too differ- 
ent from the observed fraction of "core-collapsed" clusters 
in our Galaxy (see Fig. 17 a). Of course a more serious 
comparison should take into account real cluster ages and 
the distribution of initial values for i r h (0) , which is rather 
uncertain. 

We are currently in the process of implementing exact 
dynamical integrations to handle binary-binary interac- 
tions more accurately in our simulations, using Fewbody, 
a new small-N integrator we have written. This integra- 
tor performs automatic classification of outcomes, auto- 
matic stability analysis of arbitrarily large hierarchies, au- 
tonomous integration termination, and stellar collisions. 

We are very grateful to Douglas Heggie, Steve McMil- 
lan, Simon Portegies Zwart, and Saul Rappaport for many 
helpful discussions. The three-body integrator used in this 
work is part of the Starlab software package developed 
by Piet Hut, Steve McMillan and Simon Portegies Zwart. 
This work was supported by NASA ATP Grant NAG5- 
12044 and NSF Grant AST-0206276. Some of our numeri- 
cal simulations were performed on parallel supercomputers 
at Boston University and NCSA under National Compu- 
tational Science Alliance Grant AST980014N. 
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Table 1 

Model parameters and results. 



Model f b t cc [or t dis ]/t rh (0) M(t cc )/M(0) M b (t cc )/M b (0) 



Plummcr 
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0.20 


Plummcr 
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0.90 


0.14 


Plummcr 
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20% 


120 
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0.28 


Plummcr 
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0.35 


King, W = 3 
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15 


0.25 


0.04 


King, Wo = 3 
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17 


0.00 


0.00 


King, W = 3 
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15 


0.00 


0.00 


King, W = 3 
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0.00 


0.00 
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0.00 


King, Wo = 7 
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0.00 


0.00 


King, W = 11 
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0.19 


King, W = 11 
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0.43 


0.05 


King, W = 11 
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0.20 


0.01 


King, W = 11 


20% 


38 
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0.00 
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Fig. 1. — Comparison between the use of direct three-body integrations (solid lines) and cross sections (dashed lines) in calculating the 
evolution of a Plummer model containing N = 3 X 10 s stars with 20% binaries initially. In both cases binary-binary interactions were 
turned off. The top panel shows the total mass in binaries. The middle panel shows the energy generated in binary-single interactions, as a 
fraction of the total initial binding energy of the cluster. The bottom panel shows (from top to bottom) the half-mass radius of single stars, 
the half-mass radius of binaries, and the core radius, in units of the initial half-mass radius. Time is given in units of the initial half-mass 
relaxation time. The agreement between the two methods is strong, although the energy production is slightly overestimated in the treatment 
based on cross sections, leading to divergent evolutions near core collapse (the model calculated with direct three-body integrations collapses 
earlier). 
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Fig. 2. — Evolution of the kinetic energy T, potential energy W and total conserved energy E, as well as the virial ratio 2T/\W\, for 
a Plummer model with N = 3 X 10 5 stars and 10 % binaries initially. The cluster remains very close to virial equilibrium throughout 
the integration, and energy conservation is maintained to within a few percent. Note that E is corrected for both the energy lost through 
evaporation, and the energy gained through binary-binary and binary— single interactions, so that E ^ T + W (except at t = 0). The true 
total energy of the cluster, T + W, increases significantly over time as a result of these interactions. Here we show the quantity that should 
be conserved, which we monitor (in addition to the virial ratio) for quality control purposes in all our calculations. 
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Fig. 3. — Evolution of an isolated Plummer model with N = 3 X 10 stars and 2% primordial binaries initially. The top panel shows the 
total cluster mass M and the total mass Mf, in binaries. The middle panel shows the energy released through binary-binary and binary— single 
interactions, in units of the initial binding energy of the cluster. The lower panel shows the core radius r c of the cluster, the half-mass radius 
r h g of single stars, and the half-mass radius r h b of binaries (solid lines). For comparison, the core radius of an equivalent Plummer model with 
2% primordial binaries but with all interactions turned off is also shown (short-dashed line). We see that even a primordial binary fraction 
as small as 2% can significantly delay core collapse, with t cc increasing by more than a factor 2 in this case. Also shown for comparison 
and testing is the core radius of an equivalent model where the binary-single interactions were computed with direct three-body integrations 
instead of cross sections (long-dashed line). Here again we note that the model based on direct integrations collapses slightly earlier. 
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Fig. 4. — Same as Fig. 3, but for a model with a 10% primordial binary fraction. Here the energy generated from binary-binary interactions 
clearly dominates that from binary-single interactions. We see that an isolated cluster with 10% binaries can be supported against collapse 
for about 70£ r h- Several cycles of gravothcrmal oscillations powered by primordial binaries are seen after the initial collapse. The oscillations 
in this case appear quasi-periodic with a period of roughly 50i r h(0). 
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Fig. 5. — Same as Fig. 3, but for a model with a 20% primordial binary fraction. Here the energy generated from binary— binary interactions 
is even more clearly dominant. The cluster is initially supported against collapse for about 125 t r h- After the first core collapse, gravothcrmal 
oscillations powered by primordial binaries continue up to ~ 10 3 t r h- By that time the total number of binaries has been reduced by a factor 
~ 10, but the primordial binary reservoir is still not exhausted. 
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Fig. 6. — Ratio of core collapse times with and without binary interactions as a function of the initial binary fraction /b for the Plummer 
models. The solid line shows a simple linear fit. 
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Fig. 7. — Evolution of the temperature profile near core collapse at t ~ 125 t r ^ for an isolated Plummer model with TV = 3 X 10 5 stars and 
20% primordial binaries initially (same model shown in Fig, 5). The temperature is defined by 3kT = ma^. A temperature inversion is clearly 
associated with re-expansion after core collapse. The profiles are truncated very near the center where the number of stars is small and the 
statistical noise becomes too large. The profiles during contraction (a), evolving from top to bottom, are separated by about 0.1t r h(0), and 
the profiles during re-expansion (b), from bottom to top, arc roughly 0.054,^(0) apart. The profiles shown by dotted lines in (b) are about 
4< r h(0) and 5t r h(0) after core collapse, indicating that the temperature inversion quickly disappears after re-expansion. 
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Fig. 8. — Evolution of a tidally truncated Wo = 7 King model with N = 3 X 10 s stars and 2% primordial binaries. Conventions are as in 
Fig. 3. Compared to an isolated Plummer model with the same number of stars and binaries, the evolution of this tidally truncated cluster 
to core collapse is only slightly faster. 



20 



FREGEAU, GURKAN, JOSHI, & RASIO 




MONTE CARLO SIMULATIONS OF GLOBULAR CLUSTER EVOLUTION. III. 



21 




Fig. 10. — Same as Fig. 8, but for a model with a 10% primordial binary fraction initially. Here complete disruption occurs before core 
collapse. Note also the absence of any core contraction initially. 



22 



FREGEAU, GURKAN, JOSHI, & RASIO 




Fig. 11. — Same as Fig. 8, but for a model with a 20% primordial binary fraction initially. Note the initial expansion of the core. Here 
also complete disruption occurs before core collapse. The apparent re-expansion of the core radius after t/t^lO) ~ 42 is a numerical artifact 
caused by the very small number of stars left in the cluster (our method of calculating r c picks up stars well outtside the true core). 
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Fig. 12. — Same as Fig. 10, but for a Wo = 11 King model (with a 10% primordial binary fraction). This initially much more centrally 
concentrated model undergoes deep core collapse as it runs out of binaries before complete disruption. Note also the significant initial 
expansion of the core needed to reach quasi-equilibrium in a few t r h(0). 



24 



FREGEAU, GURKAN, JOSHI, & RASIO 




Fig. 13. — Same as Fig. 10, but for a Wo = 3 King model (with a 10% primordial binary fraction). This cluster is much less centrally 
concentrated initially and therefore, as expected, it is disrupted before undergoing deep core collapse. Note, however, that significant core 
contraction occurs throughout the evolution. 
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Fig. 14. — Evolution of the binary fraction /(, in different regions for various King models. The thin solid line is for a Wo = 11 model with 
10% binaries and the thin dotted line is for a Wo = 3 model with 10% binaries. The other lines are for Wo = 7 models with increasing binary 
fractions from bottom to top as in Figs. 8—11. Binary fractions as high as 0.5—0.8 (but more typically ~ 0.1 — 0.2) can be expected in cluster 
cores, while in the outer halo one has typically < 0.1. 
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Fig. 15. — Evolution of the binary hardness and radial distributions for a Wo = 7 King model with 20% binaries. Binaries undergo clear 
mass segregation, and harden on average by about two orders of magnitude before cluster disruption. 
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Fig. 16. — Ratio of half-mass to core radius and concentration parameter c = log(rt/r c ) for various King models. The thin solid line is for 
a Wo = 11 model with 10% binaries and the thin dotted line is for a Wo = 3 model with 10% binaries. The thick lines are for Wo = 7 models 
with increasing binary fractions from top to bottom, going from 2%, 5%, 10%, to 20%. 
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Fig. 17. — Observed distribution of r^/r c and concentration parameter c for Galactic globular clusters. Clusters classified obscrvationally 
as "core-collapsed" are included in (a) and (b), but not in (c). The observed values of these basic structural parameters for the non- "core- 
collapsed" clusters are clearly in general agreement with values predicted by our simple models for clusters supported by primordial binary 
burning. 



